﻿using System;
using LatoolNet;

namespace CSExample {
  class example {

    static void Main(string[] args) {

      ComplexExample();
      MatrixInverseExample();
      MatrixFactorizeSolveExample();
      MatrixSolveExample();
      ComplexMatrixInverseExample();
      TridiagonalMatrixExample();
      PolynomialApproximationExample();
      SingularValueDecompositionExample();
      EigenproblemExample();


    }

    private static void ComplexExample() {

      Complex c1 = new Complex(1.0, 1.0);

      Complex c2 = new Complex();
      c2.Real = 2.0;
      c2.Imag = 2.0;

      Complex c3 = c1 + c2;

      // e ^ (pi * i) should give -1
      Complex c4 = Complex.Pow(Math.E, Math.PI * new Complex(0, 1));

      Console.WriteLine(c3.ToString(4));
      Console.WriteLine(c4.ToString(4));

      //Expected output
      //(3.0000,3.0000)
      //(-1.0000,0.0000)

    }

    private static void MatrixInverseExample() {

      Matrix a = new Matrix(3, 3);

      a[0, 0] = 1.0;
      a[0, 1] = -2.0;
      a[0, 2] = 3.0;

      a[1, 0] = 3.0;
      a[1, 1] = -1.0;
      a[1, 2] = 2.0;

      a[2, 0] = -1.0;
      a[2, 1] = -2.0;
      a[2, 2] = 3.0;

      Matrix ainv = a.Clone().Inv();

      Matrix unit = a * ainv;

      Console.WriteLine(unit.ToString(4));

      //Expected output
      //[1.0000,0.0000,0.0000]
      //[0.0000,1.0000,0.0000]
      //[0.0000,0.0000,1.0000]

    }

    private static void MatrixFactorizeSolveExample() {

      Matrix a = new Matrix(3, 3);

      a[0, 0] = 2;
      a[0, 1] = 3;
      a[0, 2] = -1;

      a[1, 0] = 4;
      a[1, 1] = 4;
      a[1, 2] = -3;

      a[2, 0] = -2;
      a[2, 1] = 3;
      a[2, 2] = -1;

      Matrix b = new Matrix(3, 1);

      b[0, 0] = 5;
      b[1, 0] = 3;
      b[2, 0] = 1;

      //LUFactorization.Factorize(a);
      LinearEquation.Factorize(a);

      //LUFactorization.Solve(a, b);
      LinearEquation.Solve(a, b);

      Console.WriteLine(b.ToString(3));

      //Expected output
      //[1.000]
      //[2.000]
      //[3.000]
    }

    private static void MatrixSolveExample() {
      Matrix a = new Matrix(3, 3);

      a[0, 0] = 2;
      a[0, 1] = 3;
      a[0, 2] = -1;

      a[1, 0] = 4;
      a[1, 1] = 4;
      a[1, 2] = -3;

      a[2, 0] = -2;
      a[2, 1] = 3;
      a[2, 2] = -1;

      Matrix b = new Matrix(3, 1);

      b[0, 0] = 5;
      b[1, 0] = 3;
      b[2, 0] = 1;

      //LUFactorization.Solve(a, b);
      LinearEquation.Solve(a, b);

      Console.WriteLine(b.ToString(3));

      //Expected output
      //[1.000]
      //[2.000]
      //[3.000]


    }

    private static void ComplexMatrixInverseExample() {
      ComplexMatrix a = new ComplexMatrix(3, 3);

      Random rgen = new Random();

      a[0, 0] = new Complex(rgen.NextDouble(), rgen.NextDouble());
      a[0, 1] = new Complex(rgen.NextDouble(), rgen.NextDouble());
      a[0, 2] = new Complex(rgen.NextDouble(), rgen.NextDouble());

      a[1, 0] = new Complex(rgen.NextDouble(), rgen.NextDouble());
      a[1, 1] = new Complex(rgen.NextDouble(), rgen.NextDouble());
      a[1, 2] = new Complex(rgen.NextDouble(), rgen.NextDouble());

      a[2, 0] = new Complex(rgen.NextDouble(), rgen.NextDouble());
      a[2, 1] = new Complex(rgen.NextDouble(), rgen.NextDouble());
      a[2, 2] = new Complex(rgen.NextDouble(), rgen.NextDouble());

      ComplexMatrix ainv = a.Clone().Inv();

      ComplexMatrix unit = a * ainv;

      Console.WriteLine(unit.ToString(3));

      //Expected output
      //[(1.000,0.000),(0.000,0.000),(0.000,0.000)]
      //[(0.000,0.000),(1.000,0.000),(0.000,0.000)]
      //[(0.000,0.000),(0.000,0.000),(1.000,0.000)]


    }

    private static void TridiagonalMatrixExample() {
      int bandwidth = 3;
      Matrix tri = new Matrix(5, 5, bandwidth);

      tri[0, 0] = 2;
      tri[0, 1] = -1;

      tri[1, 0] = -1;
      tri[1, 1] = 2;
      tri[1, 2] = -1;

      tri[2, 1] = -1;
      tri[2, 2] = 2;
      tri[2, 3] = -1;

      tri[3, 2] = -1;
      tri[3, 3] = 2;
      tri[3, 4] = -1;

      tri[4, 3] = -1;
      tri[4, 4] = 2;

      Console.WriteLine(tri.ToString(3));

      //Expected output
      //[2.000,-1.000,0.000,0.000,0.000]
      //[-1.000,2.000,-1.000,0.000,0.000]
      //[0.000,-1.000,2.000,-1.000,0.000]
      //[0.000,0.000,-1.000,2.000,-1.000]
      //[0.000,0.000,0.000,-1.000,2.000]

      Matrix expected_x = new Matrix(5, 1);
      expected_x[0, 0] = 1;
      expected_x[1, 0] = 3;
      expected_x[2, 0] = 3;
      expected_x[3, 0] = 3;
      expected_x[4, 0] = 1;

      Matrix b = tri * expected_x;

      //LUFactorization.Solve(tri, b);
      LinearEquation.Solve(tri, b);

      Console.WriteLine(b.ToString(3));

      //Expected output
      //[1.000]
      //[3.000]
      //[3.000]
      //[3.000]
      //[1.000]

    }

    private static void PolynomialApproximationExample() {
      double a = 1.234;
      double b = 2.354;
      double c = -4.245;
      double d = 2.987;
      //y = a * x^3 + b * x^2 + c * x + d

      double[] x = new double[5];
      double[] y = new double[5];

      x[0] = -6;
      y[0] = a * x[0] * x[0] * x[0] + b * x[0] * x[0] + c * x[0] + d;

      x[1] = -3;
      y[1] = a * x[1] * x[1] * x[1] + b * x[1] * x[1] + c * x[1] + d;

      x[2] = 0;
      y[2] = a * x[2] * x[2] * x[2] + b * x[2] * x[2] + c * x[2] + d;

      x[3] = 3;
      y[3] = a * x[3] * x[3] * x[3] + b * x[3] * x[3] + c * x[3] + d;

      x[4] = 6;
      y[4] = a * x[4] * x[4] * x[4] + b * x[4] * x[4] + c * x[4] + d;

      int dimension = 3;
      Polynomial poly = PolynomialFitting.Fit(x, y, dimension);

      double[] coefs = poly.Coeffcients;

      for (int i = coefs.Length - 1; i >= 0; i--) {
        Console.WriteLine(coefs[i].ToString("#.###"));
      }
      Console.WriteLine();
      //Expected output
      //1.234
      //2.354
      //-4.245
      //2.987
    }

    private static void SingularValueDecompositionExample() {

      Matrix mat = new Matrix(6, 4);

      mat[0, 0] = 2.27;
      mat[0, 1] = -1.54;
      mat[0, 2] = 1.15;
      mat[0, 3] = -1.94;

      mat[1, 0] = 0.28;
      mat[1, 1] = -1.67;
      mat[1, 2] = 0.94;
      mat[1, 3] = -0.78;

      mat[2, 0] = -0.48;
      mat[2, 1] = -3.09;
      mat[2, 2] = 0.99;
      mat[2, 3] = -0.21;

      mat[3, 0] = 1.07;
      mat[3, 1] = 1.22;
      mat[3, 2] = 0.79;
      mat[3, 3] = 0.63;

      mat[4, 0] = -2.35;
      mat[4, 1] = 2.93;
      mat[4, 2] = -1.45;
      mat[4, 3] = 2.30;

      mat[5, 0] = 0.62;
      mat[5, 1] = -7.39;
      mat[5, 2] = 1.03;
      mat[5, 3] = -2.57;

      Console.WriteLine(mat.ToString());

      Matrix sigma;
      Matrix U;
      Matrix VT;

      SingularValueDecomposition.Decompose(mat, out sigma, out  U, out VT);

      Console.WriteLine(sigma.ToString());
      Console.WriteLine(U.ToString());
      Console.WriteLine(VT.ToString());

      //The matrix below should be equal to the original.
      Matrix verify = U * sigma * VT;

      Console.WriteLine(verify.ToString());
    }

    private static void EigenproblemExample() {

      int rownum = 4;
      int colnum = 4;

      //Only symmetric type is supported for eigen problem. 
      Matrix mat = new Matrix(rownum, colnum, MatrixType.DoubleSymmetric);

      mat[0, 0] = 1.0;
      mat[0, 1] = 2.0;
      mat[0, 2] = 3.0;
      mat[0, 3] = 4.0;

      mat[1, 1] = 2.0;
      mat[1, 2] = 3.0;
      mat[1, 3] = 4.0;

      mat[2, 2] = 3.0;
      mat[2, 3] = 4.0;

      mat[3, 3] = 4.0;

      Matrix orig = mat.Clone();

      double[] resultValues;
      Matrix resultVectors;

      Eigenproblem.Solve(mat, out resultValues, out resultVectors);

      for (int i = 0; i < resultValues.Length; i++) {
        double lambda = resultValues[i];
        Matrix vector = Matrix.ColVector(resultVectors, i);

        //Since A x = lambda x, lhs and rhs should be same.
        Matrix lhs = orig * vector;
        Matrix rhs = lambda * vector;

        Console.WriteLine(lhs.ToString());
        Console.WriteLine(rhs.ToString());

      }

    }

  }

}
